! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Copyright (c) 2015, Regents of the University of Colorado
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are 
! permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this list of 
!    conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice, this list
!    of conditions and the following disclaimer in the documentation and/or other 
!    materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors may be 
!    used to endorse or promote products derived from this software without specific prior
!    written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 
! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 
! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL 
! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT 
! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! History:
! July 2006: John Haynes      - Initial version
! May 2015:  Dustin Swales    - Modified for COSPv2.0
! 
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
module optics_lib
  USE COSP_KINDS,     ONLY: wp
  use mod_cosp_error, ONLY: errorMessage
  implicit none

contains

  ! ##############################################################################
  !                           subroutine M_WAT
  ! ##############################################################################
  subroutine m_wat(freq, tk, n_r, n_i)
    ! ############################################################################
    !  
    ! Purpose:
    !   compute complex index of refraction of liquid water
    !
    ! Inputs:
    !   [freq]    frequency (GHz)
    !   [tk]       temperature (K)
    !
    ! Outputs:
    !   [n_r]     real part index of refraction
    !   [n_i]     imaginary part index of refraction
    !
    ! Reference:
    !   Based on the work of Ray (1972)
    !
    ! Coded:
    !   03/22/05  John Haynes (haynes@atmos.colostate.edu)
    ! ############################################################################

    ! INPUTS
    real(wp), intent(in) :: &
         freq, & ! Frequency (GHz)
         tk      ! Temperature (K)
  
    ! OUTPUTS
    real(wp), intent(out) :: &
         n_r,  & ! Real part of index of refraction
         n_i     ! Imaginary part of index of refraction

    ! Internal variables
    real(wp) :: ld,es,ei,a,ls,sg,tm1,cos1,sin1,e_r,e_i,pi,tc
    complex(wp) :: e_comp, sq

    tc = tk - 273.15_wp

    ld = 100._wp*2.99792458E8_wp/(freq*1E9_wp)
    es = 78.54_wp*(1-(4.579E-3_wp*(tc-25._wp)+1.19E-5_wp*(tc-25._wp)**2 &
         -2.8E-8_wp*(tc-25._wp)**3))
    ei = 5.27137_wp+0.021647_wp*tc-0.00131198_wp*tc**2
    a = -(16.8129_wp/(tc+273._wp))+0.0609265_wp
    ls = 0.00033836_wp*exp(2513.98_wp/(tc+273._wp))
    sg = 12.5664E8_wp
    
    tm1 = (ls/ld)**(1-a)
    pi = acos(-1._wp)
    cos1 = cos(0.5_wp*a*pi)
    sin1 = sin(0.5_wp*a*pi)
    
    e_r = ei + (((es-ei)*(1.+tm1*sin1))/(1._wp+2*tm1*sin1+tm1**2))
    e_i = (((es-ei)*tm1*cos1)/(1._wp+2*tm1*sin1+tm1**2)) &
         +((sg*ld)/1.885E11_wp)
    
!ds    e_comp = cmplx(e_r,e_i,Kind=Kind(0d0))
    e_comp = cmplx(e_r,e_i,Kind=wp)
    sq = sqrt(e_comp)
    
    n_r = real(sq)
    n_i = aimag(sq)      
    
    return
  end subroutine m_wat

  ! ############################################################################
  !                           subroutine M_ICE
  ! ############################################################################
  subroutine m_ice(freq,t,n_r,n_i)
    ! ##########################################################################
    !
    ! Purpose:
    !   compute complex index of refraction of ice
    !
    ! Inputs:
    !   [freq]    frequency (GHz)
    !   [t]       temperature (K)
    !
    ! Outputs:
    !   [n_r]     real part index of refraction
    !   [n_i]     imaginary part index of refraction
    !
    ! Reference:
    !    Fortran 90 port from IDL of REFICE by Stephen G. Warren
    !
    ! Modified:
    !   05/31/05  John Haynes (haynes@atmos.colostate.edu)
    ! ##########################################################################

    ! INPUTS
    real(wp), intent(in) :: &
         freq, & ! Frequency (GHz)
         t       ! Temperature (K)
  
    ! OUTPUTS
    real(wp), intent(out) :: &
         n_r,  & ! Real part of index of refraction
         n_i     ! Imaginary part of index of refraction

    ! Internal variables
    integer  :: i,lt1,lt2
    real(wp) :: alam,pi,t1,t2, &
         x,x1,x2,y,y1,y2,ylo,yhi,tk


    ! Parameters:
    integer,parameter :: &
         nwl  = 468,      & !
         nwlt = 62          !
    real(wp),parameter,dimension(4) :: &
         temref = [272.16,268.16,253.16,213.16]
    real(wp),parameter :: & !
         wlmin  = 0.045,  & !
         wlmax  = 8.6e6,  & !
         cutice = 167.0
    real(wp),parameter,dimension(nwlt) :: &
         wlt = &
         [0.1670e+03, 0.1778e+03, 0.1884e+03, 0.1995e+03, 0.2113e+03, 0.2239e+03, &
          0.2371e+03, 0.2512e+03, 0.2661e+03, 0.2818e+03, 0.2985e+03, 0.3162e+03, &
          0.3548e+03, 0.3981e+03, 0.4467e+03, 0.5012e+03, 0.5623e+03, 0.6310e+03, &
          0.7943e+03, 0.1000e+04, 0.1259e+04, 0.2500e+04, 0.5000e+04, 0.1000e+05, &
          0.2000e+05, 0.3200e+05, 0.3500e+05, 0.4000e+05, 0.4500e+05, 0.5000e+05, &
          0.6000e+05, 0.7000e+05, 0.9000e+05, 0.1110e+06, 0.1200e+06, 0.1300e+06, &
          0.1400e+06, 0.1500e+06, 0.1600e+06, 0.1700e+06, 0.1800e+06, 0.2000e+06, &
          0.2500e+06, 0.2900e+06, 0.3200e+06, 0.3500e+06, 0.3800e+06, 0.4000e+06, &
          0.4500e+06, 0.5000e+06, 0.6000e+06, 0.6400e+06, 0.6800e+06, 0.7200e+06, &
          0.7600e+06, 0.8000e+06, 0.8400e+06, 0.9000e+06, 0.1000e+07, 0.2000e+07, &
          0.5000e+07,0.8600e+07]
    real(wp),parameter,dimension(nwl) :: &
         tabim = &
         [0.1640e+00, 0.1730e+00, 0.1830e+00, 0.1950e+00, 0.2080e+00, 0.2230e+00, &
          0.2400e+00, 0.2500e+00, 0.2590e+00, 0.2680e+00, 0.2790e+00, 0.2970e+00, &
          0.3190e+00, 0.3400e+00, 0.3660e+00, 0.3920e+00, 0.4160e+00, 0.4400e+00, &
          0.4640e+00, 0.4920e+00, 0.5170e+00, 0.5280e+00, 0.5330e+00, 0.5340e+00, &
          0.5310e+00, 0.5240e+00, 0.5100e+00, 0.5000e+00, 0.4990e+00, 0.4680e+00, &
          0.3800e+00, 0.3600e+00, 0.3390e+00, 0.3180e+00, 0.2910e+00, 0.2510e+00, &
          0.2440e+00, 0.2390e+00, 0.2390e+00, 0.2440e+00, 0.2470e+00, 0.2240e+00, &
          0.1950e+00, 0.1740e+00, 0.1720e+00, 0.1800e+00, 0.1940e+00, 0.2130e+00, &
          0.2430e+00, 0.2710e+00, 0.2890e+00, 0.3340e+00, 0.3440e+00, 0.3820e+00, &
          0.4010e+00, 0.4065e+00, 0.4050e+00, 0.3890e+00, 0.3770e+00, 0.3450e+00, &
          0.3320e+00, 0.3150e+00, 0.2980e+00, 0.2740e+00, 0.2280e+00, 0.1980e+00, &
          0.1720e+00, 0.1560e+00, 0.1100e+00, 0.8300e-01, 0.5800e-01, 0.2200e-01, &
          0.1000e-01, 0.3000e-02, 0.1000e-02, 0.3000e-03, 0.1000e-03, 0.3000e-04, &
          0.1000e-04, 0.3000e-05, 0.1000e-05, 0.7000e-06, 0.4000e-06, 0.2000e-06, &
          0.1000e-06, 0.6377e-07, 0.3750e-07, 0.2800e-07, 0.2400e-07, 0.2200e-07, &
          0.1900e-07, 0.1750e-07, 0.1640e-07, 0.1590e-07, 0.1325e-07, 0.8623e-08, &
          0.5504e-08, 0.3765e-08, 0.2710e-08, 0.2510e-08, 0.2260e-08, 0.2080e-08, &
          0.1910e-08, 0.1540e-08, 0.1530e-08, 0.1550e-08, 0.1640e-08, 0.1780e-08, &
          0.1910e-08, 0.2140e-08, 0.2260e-08, 0.2540e-08, 0.2930e-08, 0.3110e-08, &
          0.3290e-08, 0.3520e-08, 0.4040e-08, 0.4880e-08, 0.5730e-08, 0.6890e-08, &
          0.8580e-08, 0.1040e-07, 0.1220e-07, 0.1430e-07, 0.1660e-07, 0.1890e-07, &
          0.2090e-07, 0.2400e-07, 0.2900e-07, 0.3440e-07, 0.4030e-07, 0.4300e-07, &
          0.4920e-07, 0.5870e-07, 0.7080e-07, 0.8580e-07, 0.1020e-06, 0.1180e-06, &
          0.1340e-06, 0.1400e-06, 0.1430e-06, 0.1450e-06, 0.1510e-06, 0.1830e-06, &
          0.2150e-06, 0.2650e-06, 0.3350e-06, 0.3920e-06, 0.4200e-06, 0.4440e-06, &
          0.4740e-06, 0.5110e-06, 0.5530e-06, 0.6020e-06, 0.7550e-06, 0.9260e-06, &
          0.1120e-05, 0.1330e-05, 0.1620e-05, 0.2000e-05, 0.2250e-05, 0.2330e-05, &
          0.2330e-05, 0.2170e-05, 0.1960e-05, 0.1810e-05, 0.1740e-05, 0.1730e-05, &
          0.1700e-05, 0.1760e-05, 0.1820e-05, 0.2040e-05, 0.2250e-05, 0.2290e-05, &
          0.3040e-05, 0.3840e-05, 0.4770e-05, 0.5760e-05, 0.6710e-05, 0.8660e-05, &
          0.1020e-04, 0.1130e-04, 0.1220e-04, 0.1290e-04, 0.1320e-04, 0.1350e-04, &
          0.1330e-04, 0.1320e-04, 0.1320e-04, 0.1310e-04, 0.1320e-04, 0.1320e-04, &
          0.1340e-04, 0.1390e-04, 0.1420e-04, 0.1480e-04, 0.1580e-04, 0.1740e-04, &
          0.1980e-04, 0.2500e-04, 0.5400e-04, 0.1040e-03, 0.2030e-03, 0.2708e-03, &
          0.3511e-03, 0.4299e-03, 0.5181e-03, 0.5855e-03, 0.5899e-03, 0.5635e-03, &
          0.5480e-03, 0.5266e-03, 0.4394e-03, 0.3701e-03, 0.3372e-03, 0.2410e-03, &
          0.1890e-03, 0.1660e-03, 0.1450e-03, 0.1280e-03, 0.1030e-03, 0.8600e-04, &
          0.8220e-04, 0.8030e-04, 0.8500e-04, 0.9900e-04, 0.1500e-03, 0.2950e-03, &
          0.4687e-03, 0.7615e-03, 0.1010e-02, 0.1313e-02, 0.1539e-02, 0.1588e-02, &
          0.1540e-02, 0.1412e-02, 0.1244e-02, 0.1068e-02, 0.8414e-03, 0.5650e-03, &
          0.4320e-03, 0.3500e-03, 0.2870e-03, 0.2210e-03, 0.2030e-03, 0.2010e-03, &
          0.2030e-03, 0.2140e-03, 0.2320e-03, 0.2890e-03, 0.3810e-03, 0.4620e-03, &
          0.5480e-03, 0.6180e-03, 0.6800e-03, 0.7300e-03, 0.7820e-03, 0.8480e-03, &
          0.9250e-03, 0.9200e-03, 0.8920e-03, 0.8700e-03, 0.8900e-03, 0.9300e-03, &
          0.1010e-02, 0.1350e-02, 0.3420e-02, 0.7920e-02, 0.2000e-01, 0.3800e-01, &
          0.5200e-01, 0.6800e-01, 0.9230e-01, 0.1270e+00, 0.1690e+00, 0.2210e+00, &
          0.2760e+00, 0.3120e+00, 0.3470e+00, 0.3880e+00, 0.4380e+00, 0.4930e+00, &
          0.5540e+00, 0.6120e+00, 0.6250e+00, 0.5930e+00, 0.5390e+00, 0.4910e+00, &
          0.4380e+00, 0.3720e+00, 0.3000e+00, 0.2380e+00, 0.1930e+00, 0.1580e+00, &
          0.1210e+00, 0.1030e+00, 0.8360e-01, 0.6680e-01, 0.5400e-01, 0.4220e-01, &
          0.3420e-01, 0.2740e-01, 0.2200e-01, 0.1860e-01, 0.1520e-01, 0.1260e-01, &
          0.1060e-01, 0.8020e-02, 0.6850e-02, 0.6600e-02, 0.6960e-02, 0.9160e-02, &
          0.1110e-01, 0.1450e-01, 0.2000e-01, 0.2300e-01, 0.2600e-01, 0.2900e-01, &
          0.2930e-01, 0.3000e-01, 0.2850e-01, 0.1730e-01, 0.1290e-01, 0.1200e-01, &
          0.1250e-01, 0.1340e-01, 0.1400e-01, 0.1750e-01, 0.2400e-01, 0.3500e-01, &
          0.3800e-01, 0.4200e-01, 0.4600e-01, 0.5200e-01, 0.5700e-01, 0.6900e-01, &
          0.7000e-01, 0.6700e-01, 0.6500e-01, 0.6400e-01, 0.6200e-01, 0.5900e-01, &
          0.5700e-01, 0.5600e-01, 0.5500e-01, 0.5700e-01, 0.5800e-01, 0.5700e-01, &
          0.5500e-01, 0.5500e-01, 0.5400e-01, 0.5200e-01, 0.5200e-01, 0.5200e-01, &
          0.5200e-01, 0.5000e-01, 0.4700e-01, 0.4300e-01, 0.3900e-01, 0.3700e-01, &
          0.3900e-01, 0.4000e-01, 0.4200e-01, 0.4400e-01, 0.4500e-01, 0.4600e-01, &
          0.4700e-01, 0.5100e-01, 0.6500e-01, 0.7500e-01, 0.8800e-01, 0.1080e+00, &
          0.1340e+00, 0.1680e+00, 0.2040e+00, 0.2480e+00, 0.2800e+00, 0.3410e+00, &
          0.3790e+00, 0.4090e+00, 0.4220e+00, 0.4220e+00, 0.4030e+00, 0.3890e+00, &
          0.3740e+00, 0.3540e+00, 0.3350e+00, 0.3150e+00, 0.2940e+00, 0.2710e+00, &
          0.2460e+00, 0.1980e+00, 0.1640e+00, 0.1520e+00, 0.1420e+00, 0.1280e+00, &
          0.1250e+00, 0.1230e+00, 0.1160e+00, 0.1070e+00, 0.7900e-01, 0.7200e-01, &
          0.7600e-01, 0.7500e-01, 0.6700e-01, 0.5500e-01, 0.4500e-01, 0.2900e-01, &
          0.2750e-01, 0.2700e-01, 0.2730e-01, 0.2890e-01, 0.3000e-01, 0.3400e-01, &
          0.5300e-01, 0.7550e-01, 0.1060e+00, 0.1350e+00, 0.1761e+00, 0.2229e+00, &
          0.2746e+00, 0.3280e+00, 0.3906e+00, 0.4642e+00, 0.5247e+00, 0.5731e+00, &
          0.6362e+00, 0.6839e+00, 0.7091e+00, 0.6790e+00, 0.6250e+00, 0.5654e+00, &
          0.5433e+00, 0.5292e+00, 0.5070e+00, 0.4883e+00, 0.4707e+00, 0.4203e+00, &
          0.3771e+00, 0.3376e+00, 0.3056e+00, 0.2835e+00, 0.3170e+00, 0.3517e+00, &
          0.3902e+00, 0.4509e+00, 0.4671e+00, 0.4779e+00, 0.4890e+00, 0.4899e+00, &
          0.4873e+00, 0.4766e+00, 0.4508e+00, 0.4193e+00, 0.3880e+00, 0.3433e+00, &
          0.3118e+00, 0.2935e+00, 0.2350e+00, 0.1981e+00, 0.1865e+00, 0.1771e+00, &
          0.1620e+00, 0.1490e+00, 0.1390e+00, 0.1200e+00, 0.9620e-01, 0.8300e-01]
    real(wp),parameter,dimension(nwl) :: &
         wl = &
         [0.4430e-01, 0.4510e-01, 0.4590e-01, 0.4680e-01, 0.4770e-01, 0.4860e-01, &
          0.4960e-01, 0.5060e-01, 0.5170e-01, 0.5280e-01, 0.5390e-01, 0.5510e-01, &
          0.5640e-01, 0.5770e-01, 0.5900e-01, 0.6050e-01, 0.6200e-01, 0.6360e-01, &
          0.6530e-01, 0.6700e-01, 0.6890e-01, 0.7080e-01, 0.7290e-01, 0.7380e-01, &
          0.7510e-01, 0.7750e-01, 0.8000e-01, 0.8270e-01, 0.8550e-01, 0.8860e-01, &
          0.9180e-01, 0.9300e-01, 0.9540e-01, 0.9920e-01, 0.1033e+00, 0.1078e+00, &
          0.1100e+00, 0.1127e+00, 0.1140e+00, 0.1181e+00, 0.1210e+00, 0.1240e+00, &
          0.1272e+00, 0.1295e+00, 0.1305e+00, 0.1319e+00, 0.1333e+00, 0.1348e+00, &
          0.1362e+00, 0.1370e+00, 0.1378e+00, 0.1387e+00, 0.1393e+00, 0.1409e+00, &
          0.1425e+00, 0.1435e+00, 0.1442e+00, 0.1450e+00, 0.1459e+00, 0.1468e+00, &
          0.1476e+00, 0.1480e+00, 0.1485e+00, 0.1494e+00, 0.1512e+00, 0.1531e+00, &
          0.1540e+00, 0.1550e+00, 0.1569e+00, 0.1580e+00, 0.1589e+00, 0.1610e+00, &
          0.1625e+00, 0.1648e+00, 0.1669e+00, 0.1692e+00, 0.1713e+00, 0.1737e+00, &
          0.1757e+00, 0.1779e+00, 0.1802e+00, 0.1809e+00, 0.1821e+00, 0.1833e+00, &
          0.1843e+00, 0.1850e+00, 0.1860e+00, 0.1870e+00, 0.1880e+00, 0.1890e+00, &
          0.1900e+00, 0.1910e+00, 0.1930e+00, 0.1950e+00, 0.2100e+00, 0.2500e+00, &
          0.3000e+00, 0.3500e+00, 0.4000e+00, 0.4100e+00, 0.4200e+00, 0.4300e+00, &
          0.4400e+00, 0.4500e+00, 0.4600e+00, 0.4700e+00, 0.4800e+00, 0.4900e+00, &
          0.5000e+00, 0.5100e+00, 0.5200e+00, 0.5300e+00, 0.5400e+00, 0.5500e+00, &
          0.5600e+00, 0.5700e+00, 0.5800e+00, 0.5900e+00, 0.6000e+00, 0.6100e+00, &
          0.6200e+00, 0.6300e+00, 0.6400e+00, 0.6500e+00, 0.6600e+00, 0.6700e+00, &
          0.6800e+00, 0.6900e+00, 0.7000e+00, 0.7100e+00, 0.7200e+00, 0.7300e+00, &
          0.7400e+00, 0.7500e+00, 0.7600e+00, 0.7700e+00, 0.7800e+00, 0.7900e+00, &
          0.8000e+00, 0.8100e+00, 0.8200e+00, 0.8300e+00, 0.8400e+00, 0.8500e+00, &
          0.8600e+00, 0.8700e+00, 0.8800e+00, 0.8900e+00, 0.9000e+00, 0.9100e+00, &
          0.9200e+00, 0.9300e+00, 0.9400e+00, 0.9500e+00, 0.9600e+00, 0.9700e+00, &
          0.9800e+00, 0.9900e+00, 0.1000e+01, 0.1010e+01, 0.1020e+01, 0.1030e+01, &
          0.1040e+01, 0.1050e+01, 0.1060e+01, 0.1070e+01, 0.1080e+01, 0.1090e+01, &
          0.1100e+01, 0.1110e+01, 0.1120e+01, 0.1130e+01, 0.1140e+01, 0.1150e+01, &
          0.1160e+01, 0.1170e+01, 0.1180e+01, 0.1190e+01, 0.1200e+01, 0.1210e+01, &
          0.1220e+01, 0.1230e+01, 0.1240e+01, 0.1250e+01, 0.1260e+01, 0.1270e+01, &
          0.1280e+01, 0.1290e+01, 0.1300e+01, 0.1310e+01, 0.1320e+01, 0.1330e+01, &
          0.1340e+01, 0.1350e+01, 0.1360e+01, 0.1370e+01, 0.1380e+01, 0.1390e+01, &
          0.1400e+01, 0.1410e+01, 0.1420e+01, 0.1430e+01, 0.1440e+01, 0.1449e+01, &
          0.1460e+01, 0.1471e+01, 0.1481e+01, 0.1493e+01, 0.1504e+01, 0.1515e+01, &
          0.1527e+01, 0.1538e+01, 0.1563e+01, 0.1587e+01, 0.1613e+01, 0.1650e+01, &
          0.1680e+01, 0.1700e+01, 0.1730e+01, 0.1760e+01, 0.1800e+01, 0.1830e+01, &
          0.1840e+01, 0.1850e+01, 0.1855e+01, 0.1860e+01, 0.1870e+01, 0.1890e+01, &
          0.1905e+01, 0.1923e+01, 0.1942e+01, 0.1961e+01, 0.1980e+01, 0.2000e+01, &
          0.2020e+01, 0.2041e+01, 0.2062e+01, 0.2083e+01, 0.2105e+01, 0.2130e+01, &
          0.2150e+01, 0.2170e+01, 0.2190e+01, 0.2220e+01, 0.2240e+01, 0.2245e+01, &
          0.2250e+01, 0.2260e+01, 0.2270e+01, 0.2290e+01, 0.2310e+01, 0.2330e+01, &
          0.2350e+01, 0.2370e+01, 0.2390e+01, 0.2410e+01, 0.2430e+01, 0.2460e+01, &
          0.2500e+01, 0.2520e+01, 0.2550e+01, 0.2565e+01, 0.2580e+01, 0.2590e+01, &
          0.2600e+01, 0.2620e+01, 0.2675e+01, 0.2725e+01, 0.2778e+01, 0.2817e+01, &
          0.2833e+01, 0.2849e+01, 0.2865e+01, 0.2882e+01, 0.2899e+01, 0.2915e+01, &
          0.2933e+01, 0.2950e+01, 0.2967e+01, 0.2985e+01, 0.3003e+01, 0.3021e+01, &
          0.3040e+01, 0.3058e+01, 0.3077e+01, 0.3096e+01, 0.3115e+01, 0.3135e+01, &
          0.3155e+01, 0.3175e+01, 0.3195e+01, 0.3215e+01, 0.3236e+01, 0.3257e+01, &
          0.3279e+01, 0.3300e+01, 0.3322e+01, 0.3345e+01, 0.3367e+01, 0.3390e+01, &
          0.3413e+01, 0.3436e+01, 0.3460e+01, 0.3484e+01, 0.3509e+01, 0.3534e+01, &
          0.3559e+01, 0.3624e+01, 0.3732e+01, 0.3775e+01, 0.3847e+01, 0.3969e+01, &
          0.4099e+01, 0.4239e+01, 0.4348e+01, 0.4387e+01, 0.4444e+01, 0.4505e+01, &
          0.4547e+01, 0.4560e+01, 0.4580e+01, 0.4719e+01, 0.4904e+01, 0.5000e+01, &
          0.5100e+01, 0.5200e+01, 0.5263e+01, 0.5400e+01, 0.5556e+01, 0.5714e+01, &
          0.5747e+01, 0.5780e+01, 0.5814e+01, 0.5848e+01, 0.5882e+01, 0.6061e+01, &
          0.6135e+01, 0.6250e+01, 0.6289e+01, 0.6329e+01, 0.6369e+01, 0.6410e+01, &
          0.6452e+01, 0.6494e+01, 0.6579e+01, 0.6667e+01, 0.6757e+01, 0.6897e+01, &
          0.7042e+01, 0.7143e+01, 0.7246e+01, 0.7353e+01, 0.7463e+01, 0.7576e+01, &
          0.7692e+01, 0.7812e+01, 0.7937e+01, 0.8065e+01, 0.8197e+01, 0.8333e+01, &
          0.8475e+01, 0.8696e+01, 0.8929e+01, 0.9091e+01, 0.9259e+01, 0.9524e+01, &
          0.9804e+01, 0.1000e+02, 0.1020e+02, 0.1031e+02, 0.1042e+02, 0.1053e+02, &
          0.1064e+02, 0.1075e+02, 0.1087e+02, 0.1100e+02, 0.1111e+02, 0.1136e+02, &
          0.1163e+02, 0.1190e+02, 0.1220e+02, 0.1250e+02, 0.1282e+02, 0.1299e+02, &
          0.1316e+02, 0.1333e+02, 0.1351e+02, 0.1370e+02, 0.1389e+02, 0.1408e+02, &
          0.1429e+02, 0.1471e+02, 0.1515e+02, 0.1538e+02, 0.1563e+02, 0.1613e+02, &
          0.1639e+02, 0.1667e+02, 0.1695e+02, 0.1724e+02, 0.1818e+02, 0.1887e+02, &
          0.1923e+02, 0.1961e+02, 0.2000e+02, 0.2041e+02, 0.2083e+02, 0.2222e+02, &
          0.2260e+02, 0.2305e+02, 0.2360e+02, 0.2460e+02, 0.2500e+02, 0.2600e+02, &
          0.2857e+02, 0.3100e+02, 0.3333e+02, 0.3448e+02, 0.3564e+02, 0.3700e+02, &
          0.3824e+02, 0.3960e+02, 0.4114e+02, 0.4276e+02, 0.4358e+02, 0.4458e+02, &
          0.4550e+02, 0.4615e+02, 0.4671e+02, 0.4736e+02, 0.4800e+02, 0.4878e+02, &
          0.5003e+02, 0.5128e+02, 0.5275e+02, 0.5350e+02, 0.5424e+02, 0.5500e+02, &
          0.5574e+02, 0.5640e+02, 0.5700e+02, 0.5746e+02, 0.5840e+02, 0.5929e+02, &
          0.6000e+02, 0.6100e+02, 0.6125e+02, 0.6250e+02, 0.6378e+02, 0.6467e+02, &
          0.6558e+02, 0.6655e+02, 0.6760e+02, 0.6900e+02, 0.7053e+02, 0.7300e+02, &
          0.7500e+02, 0.7629e+02, 0.8000e+02, 0.8297e+02, 0.8500e+02, 0.8680e+02, &
          0.9080e+02, 0.9517e+02, 0.1000e+03, 0.1200e+03, 0.1500e+03, 0.1670e+03]
    real(wp),parameter,dimension(nwlt,4) :: &
         tabimt = reshape(source= &
        (/0.8300e-01, 0.6900e-01, 0.5700e-01, 0.4560e-01, 0.3790e-01, 0.3140e-01, &
          0.2620e-01, 0.2240e-01, 0.1960e-01, 0.1760e-01, 0.1665e-01, 0.1620e-01, &
          0.1550e-01, 0.1470e-01, 0.1390e-01, 0.1320e-01, 0.1250e-01, 0.1180e-01, & 
          0.1060e-01, 0.9540e-02, 0.8560e-02, 0.6210e-02, 0.4490e-02, 0.3240e-02, &
          0.2340e-02, 0.1880e-02, 0.1740e-02, 0.1500e-02, 0.1320e-02, 0.1160e-02, &
          0.8800e-03, 0.6950e-03, 0.4640e-03, 0.3400e-03, 0.3110e-03, 0.2940e-03, &
          0.2790e-03, 0.2700e-03, 0.2640e-03, 0.2580e-03, 0.2520e-03, 0.2490e-03, &
          0.2540e-03, 0.2640e-03, 0.2740e-03, 0.2890e-03, 0.3050e-03, 0.3150e-03, &
          0.3460e-03, 0.3820e-03, 0.4620e-03, 0.5000e-03, 0.5500e-03, 0.5950e-03, &
          0.6470e-03, 0.6920e-03, 0.7420e-03, 0.8200e-03, 0.9700e-03, 0.1950e-02, &
          0.5780e-02, 0.9700e-02, 0.8300e-01, 0.6900e-01, 0.5700e-01, 0.4560e-01, &
          0.3790e-01, 0.3140e-01, 0.2620e-01, 0.2240e-01, 0.1960e-01, 0.1760e-01, &
          0.1665e-01, 0.1600e-01, 0.1500e-01, 0.1400e-01, 0.1310e-01, 0.1230e-01, &
          0.1150e-01, 0.1080e-01, 0.9460e-02, 0.8290e-02, 0.7270e-02, 0.4910e-02, &
          0.3300e-02, 0.2220e-02, 0.1490e-02, 0.1140e-02, 0.1060e-02, 0.9480e-03, &
          0.8500e-03, 0.7660e-03, 0.6300e-03, 0.5200e-03, 0.3840e-03, 0.2960e-03, &
          0.2700e-03, 0.2520e-03, 0.2440e-03, 0.2360e-03, 0.2300e-03, 0.2280e-03, &
          0.2250e-03, 0.2200e-03, 0.2160e-03, 0.2170e-03, 0.2200e-03, 0.2250e-03, &
          0.2320e-03, 0.2390e-03, 0.2600e-03, 0.2860e-03, 0.3560e-03, 0.3830e-03, &
          0.4150e-03, 0.4450e-03, 0.4760e-03, 0.5080e-03, 0.5400e-03, 0.5860e-03, &
          0.6780e-03, 0.1280e-02, 0.3550e-02, 0.5600e-02, 0.8300e-01, 0.6900e-01, &
          0.5700e-01, 0.4560e-01, 0.3790e-01, 0.3140e-01, 0.2620e-01, 0.2190e-01, &
          0.1880e-01, 0.1660e-01, 0.1540e-01, 0.1470e-01, 0.1350e-01, 0.1250e-01, &
          0.1150e-01, 0.1060e-01, 0.9770e-02, 0.9010e-02, 0.7660e-02, 0.6520e-02, &
          0.5540e-02, 0.3420e-02, 0.2100e-02, 0.1290e-02, 0.7930e-03, 0.5700e-03, &
          0.5350e-03, 0.4820e-03, 0.4380e-03, 0.4080e-03, 0.3500e-03, 0.3200e-03, &
          0.2550e-03, 0.2120e-03, 0.2000e-03, 0.1860e-03, 0.1750e-03, 0.1660e-03, &
          0.1560e-03, 0.1490e-03, 0.1440e-03, 0.1350e-03, 0.1210e-03, 0.1160e-03, &
          0.1160e-03, 0.1170e-03, 0.1200e-03, 0.1230e-03, 0.1320e-03, 0.1440e-03, &
          0.1680e-03, 0.1800e-03, 0.1900e-03, 0.2090e-03, 0.2160e-03, 0.2290e-03, &
          0.2400e-03, 0.2600e-03, 0.2920e-03, 0.6100e-03, 0.1020e-02, 0.1810e-02, &
          0.8300e-01, 0.6900e-01, 0.5700e-01, 0.4450e-01, 0.3550e-01, 0.2910e-01, &
          0.2440e-01, 0.1970e-01, 0.1670e-01, 0.1400e-01, 0.1235e-01, 0.1080e-01, &
          0.8900e-02, 0.7340e-02, 0.6400e-02, 0.5600e-02, 0.5000e-02, 0.4520e-02, &
          0.3680e-02, 0.2990e-02, 0.2490e-02, 0.1550e-02, 0.9610e-03, 0.5950e-03, &
          0.3690e-03, 0.2670e-03, 0.2510e-03, 0.2290e-03, 0.2110e-03, 0.1960e-03, &
          0.1730e-03, 0.1550e-03, 0.1310e-03, 0.1130e-03, 0.1060e-03, 0.9900e-04, &
          0.9300e-04, 0.8730e-04, 0.8300e-04, 0.7870e-04, 0.7500e-04, 0.6830e-04, &
          0.5600e-04, 0.4960e-04, 0.4550e-04, 0.4210e-04, 0.3910e-04, 0.3760e-04, &
          0.3400e-04, 0.3100e-04, 0.2640e-04, 0.2510e-04, 0.2430e-04, 0.2390e-04, &
          0.2370e-04, 0.2380e-04, 0.2400e-04, 0.2460e-04, 0.2660e-04, 0.4450e-04, &
          0.8700e-04, 0.1320e-03/),shape=(/nwlt,4/))

    real(wp),parameter,dimension(nwl) :: &
         tabre = &
         [0.83441,   0.83676,   0.83729,   0.83771,   0.83827,   0.84038, &
          0.84719,   0.85522,   0.86047,   0.86248,   0.86157,   0.86093, &
          0.86419,   0.86916,   0.87764,   0.89296,   0.91041,   0.93089, &
          0.95373,   0.98188,   1.02334,   1.06735,   1.11197,   1.13134, &
          1.15747,   1.20045,   1.23840,   1.27325,   1.32157,   1.38958, &
          1.41644,   1.40906,   1.40063,   1.40169,   1.40934,   1.40221, &
          1.39240,   1.38424,   1.38075,   1.38186,   1.39634,   1.40918, &
          1.40256,   1.38013,   1.36303,   1.34144,   1.32377,   1.30605, &
          1.29054,   1.28890,   1.28931,   1.30190,   1.32025,   1.36302, &
          1.41872,   1.45834,   1.49028,   1.52128,   1.55376,   1.57782, &
          1.59636,   1.60652,   1.61172,   1.61919,   1.62522,   1.63404, &
          1.63689,   1.63833,   1.63720,   1.63233,   1.62222,   1.58269, &
          1.55635,   1.52453,   1.50320,   1.48498,   1.47226,   1.45991, &
          1.45115,   1.44272,   1.43498,   1.43280,   1.42924,   1.42602, &
          1.42323,   1.42143,   1.41897,   1.41660,   1.41434,   1.41216, &
          1.41006,   1.40805,   1.40423,   1.40067,   1.38004,   1.35085, &
          1.33394,   1.32492,   1.31940,   1.31854,   1.31775,   1.31702, &
          1.31633,   1.31569,   1.31509,   1.31452,   1.31399,   1.31349, &
          1.31302,   1.31257,   1.31215,   1.31175,   1.31136,   1.31099, &
          1.31064,   1.31031,   1.30999,   1.30968,   1.30938,   1.30909, &
          1.30882,   1.30855,   1.30829,   1.30804,   1.30780,   1.30756, &
          1.30733,   1.30710,   1.30688,   1.30667,   1.30646,   1.30625, &
          1.30605,   1.30585,   1.30566,   1.30547,   1.30528,   1.30509, &
          1.30491,   1.30473,   1.30455,   1.30437,   1.30419,   1.30402, &
          1.30385,   1.30367,   1.30350,   1.30333,   1.30316,   1.30299, &
          1.30283,   1.30266,   1.30249,   1.30232,   1.30216,   1.30199, &
          1.30182,   1.30166,   1.30149,   1.30132,   1.30116,   1.30099, &
          1.30082,   1.30065,   1.30048,   1.30031,   1.30014,   1.29997, &
          1.29979,   1.29962,   1.29945,   1.29927,   1.29909,   1.29891, &
          1.29873,   1.29855,   1.29837,   1.29818,   1.29800,   1.29781, &
          1.29762,   1.29743,   1.29724,   1.29705,   1.29686,   1.29666, &
          1.29646,   1.29626,   1.29605,   1.29584,   1.29563,   1.29542, &
          1.29521,   1.29499,   1.29476,   1.29453,   1.29430,   1.29406, &
          1.29381,   1.29355,   1.29327,   1.29299,   1.29272,   1.29252, &
          1.29228,   1.29205,   1.29186,   1.29167,   1.29150,   1.29130, &
          1.29106,   1.29083,   1.29025,   1.28962,   1.28891,   1.28784, &
          1.28689,   1.28623,   1.28521,   1.28413,   1.28261,   1.28137, &
          1.28093,   1.28047,   1.28022,   1.27998,   1.27948,   1.27849, &
          1.27774,   1.27691,   1.27610,   1.27535,   1.27471,   1.27404, &
          1.27329,   1.27240,   1.27139,   1.27029,   1.26901,   1.26736, &
          1.26591,   1.26441,   1.26284,   1.26036,   1.25860,   1.25815, &
          1.25768,   1.25675,   1.25579,   1.25383,   1.25179,   1.24967, &
          1.24745,   1.24512,   1.24266,   1.24004,   1.23725,   1.23270, &
          1.22583,   1.22198,   1.21548,   1.21184,   1.20790,   1.20507, &
          1.20209,   1.19566,   1.17411,   1.14734,   1.10766,   1.06739, &
          1.04762,   1.02650,   1.00357,   0.98197,   0.96503,   0.95962, &
          0.97269,   0.99172,   1.00668,   1.02186,   1.04270,   1.07597, &
          1.12954,   1.21267,   1.32509,   1.42599,   1.49656,   1.55095, &
          1.59988,   1.63631,   1.65024,   1.64278,   1.62691,   1.61284, &
          1.59245,   1.57329,   1.55770,   1.54129,   1.52654,   1.51139, &
          1.49725,   1.48453,   1.47209,   1.46125,   1.45132,   1.44215, &
          1.43366,   1.41553,   1.39417,   1.38732,   1.37735,   1.36448, &
          1.35414,   1.34456,   1.33882,   1.33807,   1.33847,   1.34053, &
          1.34287,   1.34418,   1.34634,   1.34422,   1.33453,   1.32897, &
          1.32333,   1.31800,   1.31432,   1.30623,   1.29722,   1.28898, &
          1.28730,   1.28603,   1.28509,   1.28535,   1.28813,   1.30156, &
          1.30901,   1.31720,   1.31893,   1.32039,   1.32201,   1.32239, &
          1.32149,   1.32036,   1.31814,   1.31705,   1.31807,   1.31953, &
          1.31933,   1.31896,   1.31909,   1.31796,   1.31631,   1.31542, &
          1.31540,   1.31552,   1.31455,   1.31193,   1.30677,   1.29934, &
          1.29253,   1.28389,   1.27401,   1.26724,   1.25990,   1.24510, &
          1.22241,   1.19913,   1.17150,   1.15528,   1.13700,   1.11808, &
          1.10134,   1.09083,   1.08734,   1.09254,   1.10654,   1.14779, &
          1.20202,   1.25825,   1.32305,   1.38574,   1.44478,   1.47170, &
          1.49619,   1.51652,   1.53328,   1.54900,   1.56276,   1.57317, &
          1.58028,   1.57918,   1.56672,   1.55869,   1.55081,   1.53807, &
          1.53296,   1.53220,   1.53340,   1.53289,   1.51705,   1.50097, &
          1.49681,   1.49928,   1.50153,   1.49856,   1.49053,   1.46070, &
          1.45182,   1.44223,   1.43158,   1.41385,   1.40676,   1.38955, &
          1.34894,   1.31039,   1.26420,   1.23656,   1.21663,   1.20233, &
          1.19640,   1.19969,   1.20860,   1.22173,   1.24166,   1.28175, &
          1.32784,   1.38657,   1.46486,   1.55323,   1.60379,   1.61877, &
          1.62963,   1.65712,   1.69810,   1.72065,   1.74865,   1.76736, &
          1.76476,   1.75011,   1.72327,   1.68490,   1.62398,   1.59596, &
          1.58514,   1.59917,   1.61405,   1.66625,   1.70663,   1.73713, &
          1.76860,   1.80343,   1.83296,   1.85682,   1.87411,   1.89110, &
          1.89918,   1.90432,   1.90329,   1.88744,   1.87499,   1.86702, &
          1.85361,   1.84250,   1.83225,   1.81914,   1.82268,   1.82961]
    real(wp),parameter,dimension(nwlt,4) :: &
         tabret = reshape( &
           source =(/1.82961,   1.83258,   1.83149, &
          1.82748,   1.82224,   1.81718,   1.81204,   1.80704,   1.80250, &
          1.79834,   1.79482,   1.79214,   1.78843,   1.78601,   1.78434, &
          1.78322,   1.78248,   1.78201,   1.78170,   1.78160,   1.78190, &
          1.78300,   1.78430,   1.78520,   1.78620,   1.78660,   1.78680, &
          1.78690,   1.78700,   1.78700,   1.78710,   1.78710,   1.78720, &
          1.78720,   1.78720,   1.78720,   1.78720,   1.78720,   1.78720, &
          1.78720,   1.78720,   1.78720,   1.78720,   1.78720,   1.78720, &
          1.78720,   1.78720,   1.78720,   1.78720,   1.78720,   1.78720, &
          1.78720,   1.78720,   1.78720,   1.78720,   1.78720,   1.78720, &
          1.78720,   1.78720,   1.78720,   1.78720,   1.78800,            &
          1.82961,   1.83258,   1.83149,   1.82748,                       &
          1.82224,   1.81718,   1.81204,   1.80704,   1.80250,   1.79834, &
          1.79482,   1.79214,   1.78843,   1.78601,   1.78434,   1.78322, &
          1.78248,   1.78201,   1.78170,   1.78160,   1.78190,   1.78300, &
          1.78430,   1.78520,   1.78610,   1.78630,   1.78640,   1.78650, &
          1.78650,   1.78650,   1.78650,   1.78650,   1.78650,   1.78650, &
          1.78650,   1.78650,   1.78650,   1.78650,   1.78650,   1.78650, &
          1.78650,   1.78650,   1.78650,   1.78650,   1.78650,   1.78650, &
          1.78650,   1.78650,   1.78650,   1.78650,   1.78650,   1.78650, &
          1.78650,   1.78650,   1.78650,   1.78650,   1.78650,   1.78650, &
          1.78650,   1.78650,   1.78650,   1.78720,                       &
          1.82961,   1.83258,   1.83149,   1.82748,   1.82224,            &
          1.81718,   1.81204,   1.80704,   1.80250,   1.79834,   1.79482, &
          1.79214,   1.78843,   1.78601,   1.78434,   1.78322,   1.78248, &
          1.78201,   1.78160,   1.78140,   1.78160,   1.78220,   1.78310, &
          1.78380,   1.78390,   1.78400,   1.78400,   1.78400,   1.78400, &
          1.78400,   1.78390,   1.78380,   1.78370,   1.78370,   1.78370, &
          1.78370,   1.78370,   1.78370,   1.78370,   1.78370,   1.78370, &
          1.78370,   1.78370,   1.78370,   1.78370,   1.78370,   1.78370, &
          1.78370,   1.78370,   1.78370,   1.78370,   1.78370,   1.78370, &
          1.78370,   1.78370,   1.78370,   1.78370,   1.78370,   1.78370, &
          1.78370,   1.78400,   1.78450,                                  &
          1.82961,   1.83258,   1.83149,   1.82748,   1.82224,   1.81718, &
          1.81204,   1.80704,   1.80250,   1.79834,   1.79482,   1.79214, &
          1.78843,   1.78601,   1.78434,   1.78322,   1.78248,   1.78201, &
          1.78150,   1.78070,   1.78010,   1.77890,   1.77790,   1.77730, &
          1.77720,   1.77720,   1.77720,   1.77720,   1.77720,   1.77720, &
          1.77720,   1.77720,   1.77720,   1.77720,   1.77720,   1.77720, &
          1.77720,   1.77720,   1.77720,   1.77720,   1.77720,   1.77720, &
          1.77720,   1.77720,   1.77720,   1.77720,   1.77720,   1.77720, &
          1.77720,   1.77720,   1.77720,   1.77720,   1.77720,   1.77720, &
          1.77720,   1.77720,   1.77720,   1.77720,   1.77720,   1.77720, &
          1.77720,   1.77800/),shape=(/nwlt,4/))

    ! #####################################################################
    ! Defines wavelength dependent complex index of refraction for ice.
    ! Allowable wavelength range extends from 0.045 microns to 8.6 meter
    ! temperature dependence only considered beyond 167 microns.
    ! 
    ! interpolation is done     n_r  vs. log(xlam)
    !                           n_r  vs.        t
    !                       log(n_i) vs. log(xlam)
    !                       log(n_i) vs.        t
    !
    ! Stephen G. Warren - 1983
    ! Dept. of Atmospheric Sciences
    ! University of Washington
    ! Seattle, Wa  98195
    !
    ! Based on
    !
    !    Warren,S.G.,1984.
    !    Optical constants of ice from the ultraviolet to the microwave.
    !    Applied Optics,23,1206-1225
    !
    ! Reference temperatures are -1.0,-5.0,-20.0, and -60.0 deg C
    ! #####################################################################

    pi  = acos(-1._wp)
    n_r = 0._wp
    n_i = 0._wp
    tk  = t
    
    ! Convert frequency to wavelength (um)
    alam=3E5_wp/freq
    if((alam < wlmin) .or. (alam > wlmax)) then
       call errorMessage('FATAL ERROR(optics/optics_lib.f90:m_ice): wavelength out of bounds')
       return
    endif
    
    if (alam < cutice) then
       ! Region from 0.045 microns to 167.0 microns - no temperature depend
       do i=2,nwl
          if(alam < wl(i)) continue
       enddo
       x1  = log(wl(i-1))
       x2  = log(wl(i))
       y1  = tabre(i-1)
       y2  = tabre(i)
       x   = log(alam)
       y   = ((x-x1)*(y2-y1)/(x2-x1))+y1
       n_r = y
       y1  = log(abs(tabim(i-1)))
       y2  = log(abs(tabim(i)))
       y   = ((x-x1)*(y2-y1)/(x2-x1))+y1
       n_i = exp(y)    
    else
       ! Region from 167.0 microns to 8.6 meters - temperature dependence
       if(tk > temref(1)) tk=temref(1)
       if(tk < temref(4)) tk=temref(4)
       do i=2,4
          if(tk.ge.temref(i)) go to 12
       enddo
12     lt1 = i
       lt2 = i-1
       do i=2,nwlt
          if(alam.le.wlt(i)) go to 14
       enddo
14     x1  = log(wlt(i-1))
       x2  = log(wlt(i))
       y1  = tabret(i-1,lt1)
       y2  = tabret(i,lt1)
       x   = log(alam)
       ylo = ((x-x1)*(y2-y1)/(x2-x1))+y1
       y1  = tabret(i-1,lt2)
       y2  = tabret(i,lt2)
       yhi = ((x-x1)*(y2-y1)/(x2-x1))+y1
       t1  = temref(lt1)
       t2  = temref(lt2)
       y   = ((tk-t1)*(yhi-ylo)/(t2-t1))+ylo
       n_r = y
       y1  = log(abs(tabimt(i-1,lt1)))
       y2  = log(abs(tabimt(i,lt1)))
       ylo = ((x-x1)*(y2-y1)/(x2-x1))+y1
       y1  = log(abs(tabimt(i-1,lt2)))
       y2  = log(abs(tabimt(i,lt2)))
       yhi = ((x-x1)*(y2-y1)/(x2-x1))+y1
       y   = ((tk-t1)*(yhi-ylo)/(t2-t1))+ylo
       n_i = exp(y)
    endif
  end subroutine m_ice

  ! ############################################################################
  ! subroutine MIEINT
  ! ############################################################################
  Subroutine MieInt(Dx, SCm, Inp, Dqv, Dqxt, Dqsc, Dbsc, Dg, Xs1, Xs2, DPh, Error)
    ! ##########################################################################
    !
    !     General purpose Mie scattering routine for single particles
    !     Author: R Grainger 1990
    !     History:
    !     G Thomas, March 2005: Added calculation of Phase function and
    !     code to ensure correct calculation of backscatter coeficient
    !     Options/Extend_Source
    !
    ! ##########################################################################
    ! INPUTS
    integer, intent(in) :: &
         Inp
    real(wp),intent(in) :: &
         Dx !
    real(wp),intent(in),dimension(Inp) :: &
         Dqv
    Complex(wp),intent(in) :: &
         SCm!

    ! OUTPUTS
    Complex(wp),intent(out),dimension(InP) :: &
         Xs1,  & !
         Xs2     !
    real(wp),intent(out) :: &
         Dqxt, & !
         Dqsc, & !
         Dg,   & !
         Dbsc    !
    real(wp),intent(out),dimension(InP) :: &
         DPh
    integer :: &
         Error   !!

    ! PARAMETERS
    Integer,parameter :: &
         Imaxx   = 12000, & !
         Itermax = 30000, & ! Must be large enough to cope with the
                            ! largest possible nmx = x * abs(scm) + 15
                            ! or nmx =  Dx + 4.05*Dx**(1./3.) + 2.0
         Imaxnp = 10000     ! Change this as required
    Real(wp),parameter :: &
         RIMax=2.5,       & ! Largest real part of refractive index
         IRIMax = -2        ! Largest imaginary part of refractive index

    ! Internal variables
    Integer :: I, NStop, NmX, N, Inp2
    Real(wp)  :: Chi,Chi0,Chi1,APsi,APsi0,APsi1,Psi,Psi0,Psi1
    Real(wp),dimension(Imaxnp) :: Pi0,Pi1,Taun
    Complex(wp) :: Ir,Cm,A,ANM1,APB,B,BNM1,AMB,Xi,Xi0,Xi1,Y
    Complex(wp),dimension(Itermax) :: D
    Complex(wp),dimension(Imaxnp) :: Sp,Sm!

    ! ACCELERATOR VARIABLES
    Integer :: Tnp1,Tnm1
    Real(wp) :: Dn, Rnx,Turbo,A2
    real(wp),dimension(Imaxnp) :: S,T
    Complex(wp) :: A1
    
    If ((Dx.Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then
       Error = 1
       Return
    EndIf
    Cm = SCm
    Ir = 1 / Cm
    Y =  Dx * Cm
    If (Dx.Lt.0.02) Then
       NStop = 2
    Else
       If (Dx.Le.8.0) Then
          NStop = Dx + 4.00*Dx**(1./3.) + 2.0
       Else
          If (Dx.Lt. 4200.0) Then
             NStop = Dx + 4.05*Dx**(1./3.) + 2.0
          Else
             NStop = Dx + 4.00*Dx**(1./3.) + 2.0
          End If
       End If
    End If
    NmX = Max(Real(NStop),Real(Abs(Y))) + 15.
    If (Nmx .gt. Itermax) then
       Error = 1
       Return
    End If
    Inp2 = Inp+1
!ds    D(NmX) = cmplx(0,0,Kind=Kind(0d0))
    D(NmX) = cmplx(0,0,Kind=wp)
    Do N = Nmx-1,1,-1
       A1 = (N+1) / Y
       D(N) = A1 - 1/(A1+D(N+1))
    End Do
    Do I =1,Inp2
       Sm(I) = cmplx(0,0,Kind=wp)
!ds       Sm(I) = cmplx(0,0,Kind=Kind(0d0))
       Sp(I) = cmplx(0,0,Kind=wp)
!ds       Sp(I) = cmplx(0,0,Kind=Kind(0d0))
       Pi0(I) = 0
       Pi1(I) = 1
    End Do
    Psi0 = Cos(Dx)
    Psi1 = Sin(Dx)
    Chi0 =-Sin(Dx)
    Chi1 = Cos(Dx)
    APsi0 = Psi0
    APsi1 = Psi1
    Xi0 = cmplx(APsi0,Chi0,Kind=wp)
!ds    Xi0 = cmplx(APsi0,Chi0,Kind=Kind(0d0))
    Xi1 = cmplx(APsi1,Chi1,Kind=wp)
!ds    Xi1 = cmplx(APsi1,Chi1,Kind=Kind(0d0))
    Dg = 0
    Dqsc = 0
    Dqxt = 0
    Tnp1 = 1
    Do N = 1,Nstop
       DN = N
       Tnp1 = Tnp1 + 2
       Tnm1 = Tnp1 - 2
       A2 = Tnp1 / (DN*(DN+1._wp))
!ds       A2 = Tnp1 / (DN*(DN+1D0))
       Turbo = (DN+1._wp) / DN
!ds       Turbo = (DN+1D0) / DN
       Rnx = DN/Dx
       Psi = Tnm1*Psi1/Dx - Psi0
!ds       Psi = Dble(Tnm1)*Psi1/Dx - Psi0
       APsi = Psi
       Chi = Tnm1*Chi1/Dx       - Chi0
       Xi = cmplx(APsi,Chi,Kind=wp)
!ds       Xi = cmplx(APsi,Chi,Kind=Kind(0d0))
       A = ((D(N)*Ir+Rnx)*APsi-APsi1) / ((D(N)*Ir+Rnx)*  Xi-  Xi1)
       B = ((D(N)*Cm+Rnx)*APsi-APsi1) / ((D(N)*Cm+Rnx)*  Xi-  Xi1)
       Dqxt = Tnp1*(A + B)+ Dqxt
!ds       Dqxt = Tnp1 *      Dble(A + B)          + Dqxt
       Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc
       If (N.Gt.1) then
          Dg = Dg + (dN*dN - 1) * (ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 *(ANM1*Conjg(BNM1)) / (dN*dN - dN)
!ds          Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN)
       End If
       Anm1 = A
       Bnm1 = B
       APB = A2 * (A + B)
       AMB = A2 * (A - B)
       Do I = 1,Inp2
          If (I.GT.Inp) Then
             S(I) = -Pi1(I)
          Else
             S(I) = Dqv(I) * Pi1(I)
          End If
          T(I) = S(I) - Pi0(I)
          Taun(I) = N*T(I) - Pi0(I)
          Sp(I) = APB * (Pi1(I) + Taun(I)) + Sp(I)
          Sm(I) = AMB * (Pi1(I) - Taun(I)) + Sm(I)
          Pi0(I) = Pi1(I)
          Pi1(I) = S(I) + T(I)*Turbo
       End Do
       Psi0 = Psi1
       Psi1 = Psi
       Apsi1 = Psi1
       Chi0 = Chi1
       Chi1 = Chi
       Xi1 = cmplx(APsi1,Chi1,Kind=wp)
!ds       Xi1 = cmplx(APsi1,Chi1,Kind=Kind(0d0))
    End Do

    If (Dg .GT.0) Dg = 2 * Dg / Dqsc
    Dqsc =  2 * Dqsc / Dx**2
    Dqxt =  2 * Dqxt / Dx**2
    Do I = 1,Inp
       Xs1(I) = (Sp(I)+Sm(I)) / 2
       Xs2(I) = (Sp(I)-Sm(I)) / 2
       Dph(I) = 2 * (Xs1(I)*Conjg(Xs1(I)) + Xs2(I)*Conjg(Xs2(I))) / (Dx**2 * Dqsc)
!ds       Dph(I) = 2 * Dble(Xs1(I)*Conjg(Xs1(I)) + Xs2(I)*Conjg(Xs2(I))) / (Dx**2 * Dqsc)
    End Do
    Dbsc = 4 * Abs(( (Sp(Inp2)+Sm(Inp2))/2 )**2) / Dx**2
    Error = 0
    Return
  End subroutine MieInt
end module optics_lib
